#include "grate_simple.h"
#include "grate_constants.h"
#include <cstring>

int tarjan_find_scc_dumb(dfa * DFA, const input_data & INPUT, char * nxt_used){
  int total = DFA->total;
  int alph_size = DFA->alph_size;
  int * mtime = new int [total + 1];
  memset(mtime, 0, sizeof(int) * (total + 1));

  int * father = new int [total + 1];

  int * stack = new int [total + 1];
  int ssize = 0;

  memset(nxt_used, 0, sizeof(char) * (total + 1));

  father[1] = 0;
  mtime[1] = 1;
  int ctime = 1, cver = 1, cedge, nver, pver;
  int scc_num = 0;
  int scc_size = 0;
  while (cver > 0){
    cedge = nxt_used[cver] & LONG_MASK;
    if (cedge == 0){ 
      nxt_used[cver] ^= 16; // is in stack
      stack[++ssize] = cver;
    }
    else{ // check the son
      pver = (*DFA)[cver][cedge-1];
      if ((pver > 0) && (father[pver] == cver) && (mtime[pver] < mtime[cver])){
        mtime[cver] = mtime[pver];
        nxt_used[cver] |= 32; // improvement
      }
    }
    if (cedge >= alph_size){
      // handle vertex away

      if ((nxt_used[cver] & 32) == 0){ // no improvement
        if (stack[ssize] == cver){
          nxt_used[cver] ^= 16; // ouf of stack
          ssize--;
        }
        else{
          scc_num++;
          if (scc_num > 1){
            delete[] mtime;
            delete[] father;
            delete[] stack;
            return -1;
          }
          do{
            pver = stack[ssize--];
            nxt_used[pver] ^= 64; // in scc
            scc_size++;
            nxt_used[pver] ^= 16; // ouf of stack
          }while (pver != cver);
        }
      }

      cver = father[cver];
      continue;
    }
    nxt_used[cver]++;
    nver = (*DFA)[cver][cedge];
    if (!nver) continue;
    if (!mtime[nver]){
      ctime++;
      mtime[nver] = ctime;
      father[nver] = cver;
      cver = nver;
      continue;
    }
    else if ((nxt_used[nver] & 16) && (mtime[nver] < mtime[cver])){
      mtime[cver] = mtime[nver];
      nxt_used[cver] |= 32; // improvement
    }
  }
  delete[] mtime;
  delete[] father;
  delete[] stack;
  return scc_size;
}

int calc_vert_incycle(dfa * DFA, const input_data & INPUT, char * nxt_used){
  int total = DFA->total;
  int alph_size = DFA->alph_size;
  int * father = new int [total + 1];
  
  memset(nxt_used, 0, sizeof(char) * (total + 1));

  father[1] = 0;
  nxt_used[1] ^= 32; // used
  int cver = 1, cedge, nver;
  int vert_incycle;
  while (cver > 0){
    cedge = nxt_used[cver] & LONG_MASK;
    if (cedge == 0){
      nxt_used[cver] ^= 16; // marked
    }
    if (cedge >= alph_size){
      nxt_used[cver] ^= 16; // unmark
      cver = father[cver];
      continue;
    }
    nxt_used[cver]++;
    nver = (*DFA)[cver][cedge];
    if (!nver) continue;
    if ((nxt_used[nver] & 32) == 0){ // new vertex
      father[nver] = cver;
      nxt_used[nver] ^= 32;
      cver = nver;
    }
    else if (nxt_used[nver] & 16){
      vert_incycle = nver;
    }
  }
  delete[] father;
  return vert_incycle;
}

void tarjan_find_scc_improved(dfa * DFA, const input_data & INPUT, const int vert_incycle, char * nxt_used){
  int total = DFA->total;
  int alph_size = DFA->alph_size;
  int * mtime = new int [total + 1];
  memset(mtime, 0, sizeof(int) * (total + 1));

  int * father = new int [total + 1];

  memset(nxt_used, 0, sizeof(char) * (total + 1));
  
  father[vert_incycle] = 0;
  mtime[vert_incycle] = 1; // Here mtime stores the 
  int ctime = 1, cver = vert_incycle, cedge, nver, pver;
  while (cver > 0){
    cedge = nxt_used[cver] & LONG_MASK;
    if (cedge == 0){
      nxt_used[cver] ^= 16; // is in stack
    }
    else{ // check the son
      pver = (*DFA)[cver][cedge-1];
      if ((pver > 0) && (father[pver] == cver) && (mtime[pver] < mtime[cver])){
        mtime[cver] = mtime[pver];
        nxt_used[cver] |= 32; // improvement
      }
    }
    if (cedge >= alph_size){
      // handle vertex away

      if ((nxt_used[cver] & 32) == 0){ // no improvement
        if (cver != vert_incycle){
          nxt_used[cver] ^= 16; // ouf of stack
        }
        /*else{
          scc_num++;
        }*/
      }

      cver = father[cver];
      continue;
    }
    nxt_used[cver]++;
    nver = (*DFA)[cver][cedge];
    if (!nver) continue;
    if (!mtime[nver]){
      ctime++;
      mtime[nver] = ctime;
      father[nver] = cver;
      cver = nver;
      continue;
    }
    else if ((nxt_used[nver] & 16) && (mtime[nver] < mtime[cver])){
      mtime[cver] = mtime[nver];
      nxt_used[cver] |= 32; // improvement
    }
  }
}

double calc_grate_simple(dfa * DFA, const input_data & INPUT, const grate_config & CONFIG){
  int total = DFA->total;

  // binary encoding here: 15-mask of number of edges done, 16-is being proceeded, 32-used, 64-vertex lie in scc
  char * nxt_used = new char [total + 1];

  int i;
  int scc_size;

  if (CONFIG.IMPROVE_TARJAN){
    int vert_incycle = calc_vert_incycle(DFA, INPUT, nxt_used);
    tarjan_find_scc_improved(DFA, INPUT, vert_incycle, nxt_used);
    scc_size = 1;
    nxt_used[vert_incycle] |= 64;
    for (i=1; i<=total; i++){
      if (nxt_used[i] & 32){
        nxt_used[i] |= 64;
        scc_size++;
      }
    }
  }
  else{
    scc_size = tarjan_find_scc_dumb(DFA, INPUT, nxt_used);
  }

  if (scc_size < 0){
    INPUT.dump("Achtung! There is more than one scc!");
    delete[] nxt_used;
    return 0;
  }

  int * scc = new int [total + 1];
  memset(scc, 0, sizeof(int) * (total + 1));
  int j = 1;
  for (i=1; i<=total; i++){
    if (nxt_used[i] & 64){
      scc[i] = -j;
      j++;
    }
  }
  delete[] nxt_used;

  int scc_edgnum;
  double grate = calc_grate_for_scc(DFA, INPUT, CONFIG, scc, scc_size, scc_edgnum);
  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump("The number of strongly-connected components: 1\n");
    INPUT.dump("The number of vertices in scc: %d\n", scc_size);
    INPUT.dump("The number of edges in scc: %d\n", scc_edgnum);   
    INPUT.dump("\n");
  }
  return grate;
}

double calc_grate_for_scc (dfa * DFA, const input_data & INPUT, const grate_config & CONFIG, int * scc, int & scc_size, int & scc_edgnum){
  int i, j;
  int total = DFA->total;
  int alph_size = DFA->alph_size;
  scc_edgnum = 0;
  scc_size = 0;

  for (i=1; i<=total; i++){
    if (scc[i] >= 0) continue;
    scc_size++;
    for (j=0; j<alph_size; j++){
      if (scc[(*DFA)[i][j]] >= 0){
        (*DFA)[-scc[i]][j] = 0;
      }
      else{
        (*DFA)[-scc[i]][j] = -scc[(*DFA)[i][j]];
        scc_edgnum++;
      }
    }
  }

  delete[] scc;
  
  if (scc_size == scc_edgnum){
    return 1.0;
  }

  bool have_edge_list = false;
  int * edge1;
  int * edge2;

  dfa * DFA_SMALL;

  if (CONFIG.GEN_EDGE_LIST==1 ||
    (CONFIG.GEN_EDGE_LIST==-1 && sizeof(int)*scc_edgnum < sizeof(double)*scc_size
        && 2*sizeof(int)*scc_edgnum < sizeof(int)*(DFA->size)*alph_size)){
    
    have_edge_list = true;
    edge1 = new int [scc_edgnum];
    edge2 = new int [scc_edgnum];
    int cnum = 0;
    for (i=1; i<=scc_size; i++){
      for (j=0; j<alph_size; j++){
        if ((*DFA)[i][j] > 0){
          edge1[cnum] = i;
          edge2[cnum++] = (*DFA)[i][j];
        }
      }
    }

    delete DFA;
  }
  else{
    DFA_SMALL = create_dfa(INPUT, scc_size + 1);
    for (i=1; i<=scc_size; i++){
      for (j=0; j<alph_size; j++){
        (*DFA_SMALL)[i][j] = (*DFA)[i][j];
      }
    }

    delete DFA;
  }

  double * cnt_old = new double [scc_size + 1];
  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump("old counters allocated\n");
  }
  double * cnt_new = new double [scc_size + 1];
  if (CONFIG.OUTPUT_PROGRESS){
    INPUT.dump("new counters allocated\n");
  }
  
  for (i=1; i<=scc_size; i++){
    cnt_new[i] = 0;
    cnt_old[i] = 1;
  }
  
  double minrate, maxrate, crate;
  int iter = 0;
  while (iter < CONFIG.MAXITER){
    iter++;
    if (have_edge_list){
      for (i=0; i<scc_edgnum; i++){
        cnt_new[edge2[i]] += cnt_old[edge1[i]];
      }
    }
    else{
      for (i=1; i<=scc_size; i++){
        for (j=0; j<alph_size; j++){
          if ((*DFA_SMALL)[i][j] > 0){
            cnt_new[(*DFA_SMALL)[i][j]] += cnt_old[i];
          }
        }
      }
    }
    minrate = 1e100;
    maxrate = 0;
    for (i=1; i<=scc_size; i++){
      if (CONFIG.ADD_LOOPS){
        cnt_new[i] += cnt_old[i];
      }
      crate = cnt_new[i] / cnt_old[i];
      if (crate < minrate) minrate = crate;
      if (crate > maxrate) maxrate = crate;
      cnt_old[i] = cnt_new[i];
      cnt_new[i] = 0;
    }
    if (maxrate - minrate < 2 * INPUT.delta){
      double res = (minrate + maxrate) / 2 - (CONFIG.ADD_LOOPS == true);
      if (minrate > MAX_GRATE) res = 0;
      if (CONFIG.OUTPUT_PROGRESS) INPUT.dump("\n%d iterations used\n\n", iter);
      if (GUI_VERSION) printf("%d\n", iter);

      if (have_edge_list){
        delete[] edge1;
        delete[] edge2;
      }
      else{
        delete DFA_SMALL;
      }
      delete[] cnt_old;
      delete[] cnt_new;

      return res;
    }
  }
  if (CONFIG.OUTPUT_PROGRESS) INPUT.dump("\nNot enough iterations or precision error! %.20lf precision only!\n\n", (maxrate-minrate)/2);
  if (GUI_VERSION) printf("-1\n");
  double res = (maxrate + minrate) / 2 - (CONFIG.ADD_LOOPS == true);

  if (have_edge_list){
    delete[] edge1;
    delete[] edge2;
  }
  else{
    delete DFA_SMALL;
  }
  delete[] cnt_old;
  delete[] cnt_new;

  return res;
}